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Abstract 

Application  of  wavelet  networks  for  identification  of  a  direct  internal  reforming  solid  oxide  fuel  cell  (DIR-SOFC)  stack  is  reported  in  this  paper. 
The  SOFC  is  a  complex  system  particularly  when  it  is  directly  fueled  with  hydrocarbons  (natural  gas,  coal  gas,  etc.).  Most  of  the  traditional  models 
of  the  SOFC,  based  on  the  reforming,  electrochemical  and  thermal  modeling,  are  too  complicated.  To  facilitate  controller  design  and  analysis  of 
systems,  the  wavelet  network  dynamic  model  of  the  DIR-SOFC  is  constructed,  avoiding  the  consideration  of  the  complex  processes  in  the  fuel  cells. 
The  input  and  output  data  are  used  for  initializing  and  training  the  wavelet  network  by  a  recursive  approach.  The  Gram-Schmidt  algorithm,  the 
Cross-Validation  method  and  immune  selection  principles  are  applied  to  optimization  of  the  network.  The  simulation  is  performed  and  comparisons 
of  characteristics  under  different  operating  conditions  are  given.  The  results  show  high  static  and  dynamic  accuracy  of  the  identified  model.  Further, 
the  obtained  wavelet  network  model  can  be  used  for  developing  the  model-based  controllers  of  DIR-SOFC. 

©  2008  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  solid  oxide  fuel  cell  (SOFC)  is  based  on  a  solid-state 
ion-conducting  electrolyte,  which  functions  at  high  tempera¬ 
ture.  Due  to  high  efficiency,  high  reliability,  and  low  levels  of 
noise  and  pollution,  the  SOFC  has  been  considered  as  one  of 
the  most  promising  technologies  for  electrical  energy  genera¬ 
tion  [1],  The  high  operating  temperature  (up  to  1000  °C)  allows 
internal  reforming  and  promotes  rapid  kinetics  with  nonprecious 
materials.  Therefore,  the  SOFC  can  be  directly  fueled  with  pure 
hydrogen,  natural  gas  and  other  hydrocarbons  [1,2], 

Due  to  a  lot  of  coupling  parameters,  the  SOFC  is  con¬ 
sidered  a  complicated  nonlinear  multi-input  and  multi-output 
(MIMO)  system  which  is  difficult  to  model.  In  the  past  cou¬ 
ple  of  years,  several  models  [3-10]  have  been  developed  and 
tested  to  study  characteristics  of  the  SOFC  stack.  Achenbach 
presented  a  mathematical  model  of  a  planar  SOFC  and  showed 
the  distributions  of  the  current  density  and  temperature  obtained 
by  computational  method  [3].  Further,  the  dynamic  voltage 
response  to  an  electrical  current  change  was  also  discussed. 
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Costamagna  and  Honegger  [4]  presented  and  validated  a  simu¬ 
lation  model  for  a  SOFC  stack  integrated  with  an  air  pre-heater. 
Recknagle  et  al.  [8]  developed  three-dimensional  thermo-fluid 
electrochemical  models  of  planar  SOFC  stacks.  With  the  aid 
of  a  simulation  tool  combining  a  computational  fluid  dynam¬ 
ics  (CFD)  code  and  an  electrochemistry  calculation  method,  the 
distributions  of  the  temperature,  current  density  and  fuel  species 
were  shown  and  investigated.  In  Refs.  [9,10],  a  one-dimensional 
model  has  been  used  for  simulation  of  an  anode- supported 
intermediate  temperature  DIR-SOFC  stack.  Further,  using  this 
model,  the  one-dimensional  steady-state  parameter  distribution 
and  dynamic  responses  to  several  current  density  step-changes 
were  generated  and  discussed.  Most  of  the  models  are  based  on 
the  fundamentals  of  heat,  momentum  and  mass  transfer,  and  are 
focused  on  explaining  the  operation  mechanism  of  the  SOFC. 
Although  the  models  describe  the  physical  and  chemical  pro¬ 
cesses  well,  most  of  them  are  too  complicated  to  be  applied  to 
controller  design. 

To  simplify  the  fuel  cell  models  for  the  performance  analysis 
and  valid  control  strategy  design,  statistical  identification  meth¬ 
ods  based  on  measurements  have  gained  attention  during  the 
recent  years  [11-13],  Arriagada  et  al.  [11]  developed  an  artificial 
neural  network  (ANN)  model  of  SOFC  that  is  applicable  to  the 
prediction  of  the  static  current  density-voltage  characteristics 
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Nomenclature 

a 

dilation  parameter  of  wavelet 

Ai,A2 

anodic  and  cathodic  reaction  surfaces 

b 

translation  parameter  of  wavelet 

B 

stoichiometric  matrix 

cP 

specific  heat  capacity  (Jkg-1  K-1) 

E 

Nemst  voltage  (V) 

E0 

EMF  (electro-motive-force)  at  standard  pressure 

£act 

activation  energy  (J  mol-1) 

F 

gas  flow  rate  (kg  s-1) 

F 

gas  flow  rate  vector 

E]  Far 

Faraday  constant  (96,485  C  mol-1) 

h 

gas  enthalpy  (Jkg-1) 

AH 

enthalpy  change  (J  mol-1) 

I 

current  density  (A  m-2) 

k 

exchange  current  density  (Am-2) 

II 

limiting  current  density  (A  m-2) 

Im 

mean  current  density  (Am-2) 

kn 

rate  coefficient  for  reforming  reaction 

Kn 

equilibrium  constant 

adsorption  constant 

Aw 

number  of  multi-dimensional  wavelets 

P 

partial  pressure  (bar) 

q 

heat  generation  (J  m-2  s-1) 

r 

reaction  rate  (molm-2  s-1) 

r 

reaction  rate  vector 

R 

gas  constant  (8.314  JK-1  mol-1) 

sb 

heat  source  (J  m-3  s-1) 

T 

temperature  (K) 

Tin 

inlet  temperature  (K) 

Tout 

outlet  temperature  (K) 

V 

gas  velocity  vector 

Vcell 

cell  voltage  (V) 

Vout 

stack  output  voltage  (V) 

UJ 

coefficient  of  wavelet 

W 

molar  mass  vector 

Greek  symbols 

am,  «c: 

1  charge  transfer  coefficients 

ym,  rCi 

1  coefficients  for  1™  and  (Am-2) 

s 

thickness  or  depth  (m) 

z 

thermal  conductivity  (W  m-1  K-1) 

p 

number  of  moles  of  electrons  participating  in  the 

reaction 

p 

density  (kg  m-3) 

f 

wavelet  function 

multi-dimensional  wavelet  function 

n 

ohmic  resistance  (Q  m2) 

Superscripts 

an 

anode 

c 

cell  unit 

ca 

cathode 

e 

electrochemical  reaction 

Is 

lower  separator 

re  reforming  reaction 

sh  water-gas  shift  reaction 

us  upper  separator 


under  different  operating  conditions.  In  Ref.  [12],  a  RBF  neural 
network  (NN)  identification  technology  was  employed  to  estab¬ 
lish  a  dynamic  model  of  a  molten  carbonate  fuel  cell  (MCFC) 
stack  that  was  used  for  designing  a  temperature  controller.  In 
Ref.  [  1 3] ,  the  Hammerstein  nonlinear  system,  composed  of  a  lin¬ 
ear  subsystem  and  a  nonlinear  subsystem,  was  applied  to  model 
the  SOFC  for  dynamic  response  studies.  The  NN  and  Ham¬ 
merstein  model  are  both  black-box  modeling  methods  which 
are  based  on  input  and  output  data,  and  need  not  the  priori 
knowledge  about  the  internal  structure  of  the  object. 

Among  the  black-box  methods,  the  wavelet  network  is  a 
novel  and  attractive  modeling  tool.  The  wavelet  theory  has  been 
applied  to  many  scientific  areas  such  as  signal  processing  and 
system  identification  [14,15].  The  wavelet  decomposition  and 
multi-resolution  approximation  (MRA)  guarantee  that  any  func¬ 
tion  of  L2(R )  can  be  approximated  to  any  prescribed  accuracy 
with  a  finite  sum  of  wavelets  [14,15].  Wavelet  networks  combine 
the  advantages  of  wavelets  and  neural  networks,  and  provide 
an  efficient  constructive  method  for  the  implementation  of  net¬ 
works  [16].  Therefore,  wavelet  networks  can  be  considered  as 
an  alternative  to  neural  and  radial  basis  function  networks.  Many 
researches  [15-19]  on  the  successful  application  of  wavelet  net¬ 
works  indicate  that  the  approach  is  interesting  and  powerful. 
Combined  with  advanced  control  methods,  wavelet  networks 
can  be  used  for  controlling  complicated  coupling  systems.  Wai 
[19]  successfully  developed  a  robust  controller  with  a  wavelet 
network  uncertainty  observer  to  control  the  slider  position  of 
a  motor-mechanism  coupling  system.  Lin  et  al.  [20]  presented 
an  adaptive  wavelet  network  control  system  to  control  a  syn¬ 
chronous  motor  servo  drive  system,  and  experimental  results 
show  an  enhanced  robust  control  performance. 

In  this  paper,  the  aim  is  to  identify  a  nonlinear  black-box 
MIMO  model  for  a  DIR-SOFC  stack  using  wavelet  networks. 
The  constructive  process  of  the  wavelet  network  mainly  consists 
of  three  parts  including  generation  of  wavelet  libraries,  evolution 
of  coefficients  and  selection  of  multi-dimensional  wavelets  used 
in  the  reconstruction  of  nonlinear  functions.  One  of  the  advan¬ 
tages  of  the  wavelet  network  is  that  it  can  be  initialized  according 
to  the  input  and  output  data  [16,17],  To  implement  the  con¬ 
struction  and  optimization  of  the  wavelet  network,  the  recursive 
approach,  the  Gram-Schmidt  algorithm,  the  Cross-Validation 
method  and  immune  selection  principles  are  employed. 

The  present  study  is  organized  into  five  sections.  In  Section 

2,  the  physical  model  of  the  DIR-SOFC  is  described.  In  Section 

3,  the  procedure  for  the  wavelet  network  modeling  is  explained. 
In  Section  4,  nonlinear  dynamic  modeling  of  the  DIR-SOFC 
stack  using  wavelet  networks  is  detailedly  described.  Further,  the 
modeling  results  and  comparison  of  static  and  dynamic  charac¬ 
teristics  under  different  operating  conditions  are  presented  and 
discussed.  Section  5  concludes  the  paper. 
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Fig.  1.  Schematic  diagram  of  a  DIR-SOFC. 


2.  Physical  model  of  DIR-SOFC 

2.1.  Reforming  model 

The  layout  of  a  DIR-SOFC  is  illustrated  in  Fig.  1.  At  the 
anode  of  the  DIR-SOFC,  the  methane  is  reformed  to  hydrogen 
for  oxidation  in  the  electrochemical  reaction  [1,2],  The  principal 
reactions  in  the  methane  steam-reforming  are 


CH4  +  2H20  4*  C02+4H2 

(D 

CH4  +  H20  4*  CO  +  3H2 

(2) 

CO  +  h2o  4*  co2  +  h2 

(3) 

The  reversible  reactions  obey  the  laws  of  chemical  equilib¬ 
rium.  We  denote  by  r\,  r2  and  7-3  the  rates  of  reactions  (1),  (2) 
and  (3),  respectively.  They  are  computed  by 

_  (M/p^)(Pch4  p^2o  -  {PAUl  PcoJK,  )) 

(4) 

n  DEN2 

(k2/P^)(PcH,PH2o  ~  (P^Pco/Ki)) 
n  ~  DEN2 

(5) 

(£3/Fh2)(FcoFh2o  -  (Ph2Pco2/k3)) 
r3  DEN2 

(6) 

DEN  is  defined  as 

DEN  1+/rC0^C0+ </JH2  +  AadH4PCH4 

o 

£ 

o 

+ 

Pn, 

(7) 

(1/2)02  +  2e-  —4  O2- 


(9) 


The  relation  between  the  cell  output  voltage  and  the  current 
density  [4,9]  is  expressed  as 


Vceii  =  E  -  /  (a  + 


aanFFar/on  aCaFFar/oa 


(10) 


where  I  is  the  current  density,  £2  the  ohmic  resistance  that 
depends  on  the  material  thickness  and  on  the  operating  tempera¬ 
ture,  R  the  gas  constant  (8.314  JK-1  mol-1),  T0  the  temperature 
of  the  cell  unit,  Fpar  the  Faraday  constant  (96,485  C  mol-1),  a3" 
and  aca  are  the  charge  transfer  coefficients  for  the  anode  and 
cathode,  respectively,  fi  is  the  number  of  moles  of  electrons  par¬ 
ticipating  in  the  reaction,  and  Ip  is  the  limiting  current  density. 
The  Nemst  voltage  of  the  cell,  E,  is  calculated  by 


RTC 

E  =  E0  +  — —  In 

2/Tar 


H2  ro2 
ph2o 


(11) 


where  £o  is  the  EMF  (electro-motive-force)  at  standard  pressure 

[1]. 

Exchange  current  densities  at  the  anode  and  cathode  are 
respectively  given  by 


where  yan  and  yca  are  coefficients,  E™t  and  E%at  are  the  activa¬ 
tion  energies,  and  Pref  is  the  reference  pressure  (1  bar)  [4]. 

The  electrochemical  reaction  rates  are  directly  related  to  the 
current  density.  The  rates  of  reactions  (8)  and  (9)  (>%  and  r$  )  are 
calculated  by 


r4  =  r5  = 


2  FFal' 


(14) 


2.3.  Thermal  model 


In  Eqs.  (4)-(7),  P  is  the  partial  pressure;  kn,  with  n  =  1-3, 
are  the  rate  coefficients;  Kn,  with  n  =  1-3,  are  the  equilibrium 
constants  for  the  reforming  reactions;  /Tad  the  adsorption  con¬ 
stants.  kn,  Kn  and  A"ad  are  the  functions  of  temperature  which 
are  presented  in  Ref.  [21]  detailedly.  Reforming  reaction  rates 
are  determined  by  the  temperature,  gas  partial  pressures,  equi¬ 
librium  constants,  catalyst,  etc.  All  these  factors  will  be  changed 
with  the  flows  and  reactions  in  the  stack. 


Reforming  reactions  are  highly  endothermic  reactions.  The 
electrochemical  reaction  is  exothermic  and  supplies  heat  and 
steam  to  the  endothermic  reforming  reactions.  Non-uniform 
reaction  rate  distribution  in  the  stack  makes  the  differences  of 
released  heat  across  the  stack,  and  the  heat  is  transferred  between 
each  component  of  the  fuel  cells  (including  convection,  radiation 
and  conduction). 

For  the  fuel  and  air  flows,  the  energy  equation  is  given  as 


2.2.  Electrochemical  model 

At  the  anode,  the  main  electrochemical  reaction  is  the  oxi¬ 
dation  of  hydrogen  [2],  At  the  cathode,  molecular  oxygen  is 
reduced  to  negatively  charged  ions. 

H2  +  O2-  H20  +  2e2-  (8) 


d(ph) 

idi-  ' ' 


+  div 


=  Sh, 


(15) 


where  p  is  the  density,  h  the  gas  enthalpy,  V  the  gas  velocity 
vector,  /  the  thermal  conductivity,  div  denotes  the  divergence, 
V  denotes  the  gradient,  Cp  the  specific  heat  capacity  and  »S'h  is 
the  heat  source. 
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For  the  cell  unit  and  the  separators,  the  energy  equation  can 
be  simplified  as 


matrices  B\  and  B2,  and  reaction  rate  vector  r  are  respectively 
given  by 


pCp—  -  div(fVr)  =  Sh.  (16) 

The  local  heat  released  by  electrochemical  reactions  is 
expressed  by 

=  -A Her4  -  Fcell/.  (17) 

where  VCeit/  is  the  local  output  electrical  power  produced  by  the 
SOFC. 

The  heat  changes  caused  by  reforming  reactions  and  the 
water-gas  shift  reaction  are  respectively  given  by 


-10  0 
0  1  0 
1  -1  0 
3  1  -1 

-1  -1  1 


r  =  (rj,  r2,  r3,  r4)T. 


(25) 


(26) 

(27) 


qle  =  —  AH\r\  —  AH2r2, 


(18)  3.  Wavelet  network  modeling 


qsh  =  -AH3r3. 


(19)  3.7.  Wavelet  decomposition 


Besides  the  heat  generation,  the  convective  and  radiant  heat 
transfers  between  the  cell  unit,  separators  and  flows  are  taken 
into  account.  For  each  object,  Sh  in  Eqs.  (15)  and  (16)  will  be 
given  differently. 

cell  unit: 

K  =  (20) 

bc 

separators: 

Sjf  =  Sjf  =  0  (21) 

fuel  flow: 


A  function  ■^r{x)eL2{R)  is  admissible  as  a  wavelet,  if  its 
Fourier  transform  i/r(&>)  satisfies 

„  f°°  mco)\2 

C,/,  =  /  - d co  <  00.  (28) 

Jo  w 

For  any  function  fix')  e  L2(R),  its  continuous  wavelet  trans¬ 
form  is  defined  as 

w(a,b)=  [  f{x)fa^b{x)dx,  (29) 

Jr 

where  aeR+  and  be  R  are  dilation  and  translation  parame¬ 
ters,  respectively  [14,15],  Jra^{x)  is  obtained  by  scaling  mother 
wavelet  fix)  by  a  and  translating  it  by  b: 


sfel 


qsh 

,5fuel 


(22) 


fa,b{x)  =  a  1/2t 


(30) 


air  flow: 

Sf  =  0  (23) 

In  Eqs.  (20)  and  (22),  Sc  and  5fuel  are  the  cell  unit  thickness 
and  the  fuel  channel  depth,  respectively. 

2.4.  Mass  conservation 

The  continuous  reforming  and  electrochemical  reactions 
change  the  chemical  species  concentrations  in  the  flows,  when 
the  fuel  and  air  pass  over  the  electrodes  of  the  SOFC.  Since 
all  gases  are  assumed  to  be  ideal  gases,  the  gas  flow  rates  are 
expressed  by 

f-rl  +  [  1 

TTj  [W2ffA2B2r5dA2\  ' 

where  Fo  is  the  inlet  gas  flow  rate,  A\  and  A3  are  anodic  and 
cathodic  reaction  surfaces,  respectively.  W\  and  W2  are  the  molar 
mass  vectors  of  the  fuel  and  air,  respectively.  Stoichiometric 


The  function  /(x)  can  be  reconstructed  by  the  inverse  wavelet 
transform 

/(x)  =  C^1  J  J  w(a,  b)ira^{xy^  db.  (31) 

As  a  matter  of  fact,  the  continuous  wavelet  transform  and  its 
inverse  transform  need  to  be  implemented  on  digital  computers. 
The  reconstruction  equation  is  thus  discretized  into 

fix)  =  (32) 

1 

where 

flix)  =  a™/2if(a™x  -  nbo),  (33) 

and  l  represents  the  integer  pair  ( m ,  n). 

For  multivariable  nonlinear  coupling  systems,  we  construct 
multi-dimensional  wavelets  by  using  the  product  of /Vjnput  (the 
number  of  input  variables)  scalar  wavelets. 
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Fig.  2.  Diagram  of  the  multi-dimensional  wavelet  reconstruction. 


where 

au  =  a  omi,  (35) 

bu  =  nia0m‘bo,  (36) 

x  =  (Xl,  (37) 

The  reconstruction  of  the  multivariable  function  can  be  writ¬ 
ten  as 

Nw 

fix)  =  Y^wMix).  (38) 


The  multi-dimensional  wavelet  decomposition  of  fix )  can  be 
regarded  as  a  one-hidden-layer  network  (shown  in  Fig.  2)  with 
'P  as  the  activation  function  of  the  hidden  neurons  and  with  the 
coefficient  wi  as  the  weight  of  the  hidden  layer. 

3.2.  Initialization  of  wavelet  networks 

Using  MRA  of  the  wavelet  transform,  it  can  be  shown  that 
the  bigger  Nw  (the  number  of  multi-dimensional  wavelets  used 
in  reconstruction)  and  wider  ranges  of  dilation  and  translation 
parameters  would  result  in  a  better  approximation  of  the  function 
[14,15],  In  practice,  however,  the  multi-dimensional  available 
data  are  finite  and  sparse.  Many  dilated  and  translated  wavelets 
do  not  contain  any  data  point  in  their  supports  and  actually  have 
no  contribution  to  the  function  reconstruction  [16].  They  will  not 
improve  the  approximation  accuracy,  but  increase  the  computa¬ 
tional  time  and  the  storage  requirements.  Therefore,  generating 
wavelet  libraries  and  eliminating  the  redundant  wavelets  are 
beneficial  and  necessary. 

3.2.1.  Generating  wavelet  libraries 

Discrete  dyadic  wavelets  (ao=2,  fro  =  1 )  are  the  typical 
biorthogonal  wavelet  bases.  We  generate  the  libraries  accord¬ 
ing  to  the  distribution  of  the  available  training  data,  and  select 
the  scalar  wavelets  from  the  libraries  to  construct  the  multi¬ 
dimensional  wavelets.  We  denote  by  [min,,  max,]  the  domain 


containing  the  values  of  the  ith  component  of  the  input  vectors. 
In  order  to  guarantee  that  the  wavelets  extend  initially  over  the 
whole  domain  of  the  ith  component  of  the  input  vectors,  relations 


(39)  and  (40)  should  be  satisfied  [17]. 

2~mi  <  0(max,  —  min,),  (39) 

and 

min,-  <  2~m’ni  <  max;,  (40) 

where  9  is  a  parameter  adjusted  according  to  the  mother  wavelet. 
Hence,  we  can  obtain 

m;  >  —  log2(fr(max;  —  min,)),  (41) 

and 

2m‘  min,-  <  «,-  <  2m  max,-.  (42) 


Under  above  two  conditions  (relations  (41)  and  (42)),  we 
generate  two  integer  sets 

Smj  =  {[— log2(fr(max,-— min,-))],  r-log2(fr(max,-min,-))]-|-l, 
f-  log2  (<9(max,-  -  min,- ))]  +2,  f-log2  (fr(max,  -min,- ))] 


+3],  (43) 

and 

Sm  =  {r2m''min,-],  [2'"!'min,-]  +  1,  ...,  [2m(max,-J},  (44) 


where  [  ]  represents  rounding  toward  ceiling,  and  [  J  represents 
rounding  toward  floor. 

Sm  and  Sm  are  the  dilation  and  translation  sets,  respectively. 
As  can  be  seen  from  Eq.  (44),  the  translations  depend  on  the 
values  of  dilations;  the  number  of  translations  increases  expo¬ 
nentially  with  nij .  We  denote  by  { Smi ,  5„;]  the  wavelet  library 
for  the  ith  component  of  the  input  vectors.  The  wavelet  library 
for  all  Mnput  input  components  is  [J {Sm;,  Sni). 

3.2.2.  Selection  of  wavelets 

From  all  the  possible  combinations  of  scalar  wavelets,  we 
select  a  number  of  combinations,  which  contain  most  useful 
information,  to  construct  the  wavelet  network.  In  Section  3.2.1, 
the  wavelet  libraries  are  generated  based  on  the  input  data,  not 
on  the  output  data.  To  eliminate  the  redundant  wavelets  in  the 
libraries  for  the  estimation  off  we  employ  the  Gram-Schmidt 
algorithm  relying  on  the  input  and  output  data. 

The  Gram-Schmidt  algorithm  is  a  recursive  process  which 
aims  to  obtain  an  orthogonal  basis  from  a  non-orthogonal 
set  [22],  It  first  selects  a  multi-dimensional  wavelet  from  the 
libraries  that  best  fits  the  training  data  as  the  initial  basis,  and  then 
repeatedly  selects  the  multi-dimensional  wavelets  in  the  remain¬ 
der  to  enlarge  the  basis  so  that  the  basis  can  best  fit  the  data.  The 
reader  is  referred  to  [22]  for  the  details  of  the  Gram-Schmidt 
algorithm. 

Assume  that  we  obtain  Nw  multi-dimensional  wavelets  by 
Gram-Schmidt  algorithm  as  follows 

[<fq,  ...<%w}. 


(45) 
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The  hidden  layer  of  the  network  is  made  up  of  the  selected 
wavelets.  The  determination  of  Nw  will  be  given  in  Section  3.3.2. 

3.3.  Training  algorithm 

3.3.1.  Coefficients  of  wavelets 

The  MIMO  system  can  be  considered  as  the  synthesis  of  sev¬ 
eral  MISO  subsystems.  Thus,  we  can  separately  train  each  MISO 
subsystem.  Considering  a  training  data  set  {(x(fc),  y(k))\k  = 
1 ,  ....  Nt]  generated  by  a  system  with  N\npm  inputs  and  a  scalar 
output,  we  use  the  recursive  method  to  update  the  coefficients. 

Iterative  updating  formulas  for  wis  (k),  which  is  the  coefficient 
of  i Pifxfk))  at  iteration  k,  are  given  as  follows: 

JVW 

e(k)  =  y(k )  -  -  mi(xfk)\  (46) 


dts(k)  =  kdifk  -  1)  +  %Cx(k)), 

Vh{x(k))e{k) 


wis(k)  =  wis(k  -  1}  +  - 


difk) 


(47) 

(48) 


In  Eq.  (47),  0  <  X  <  1  is  the  forgetting  factor  that  indicates 
the  different  importance  of  data.  The  data  near  iteration  k  are 
more  important  than  those  far  away  [23].  The  derivation  of  Eqs. 
(46)-(48)  is  presented  in  Appendix  A.  Updating  recursions  for 
coefficients  of  other  wavelets  are  similar  to  Eqs.  (46)-(48). 

3.3.2.  Number  of  multi-dimensional  wavelets 

In  Section  3.2.2,  we  assume  the  number  of  multi-dimensional 
wavelets  in  the  network  is  Nw,  and  select  the  wavelets  by 
Gram-Schmidt  algorithm.  However,  choosing  the  number  Nw 
is  a  difficult  process.  The  wavelet  network  model  is  expected 
to  have  relatively  high  accuracy  when  evaluated  not  only  with 
training  data  but  also  with  fresh  data.  If  Nw  is  over  large,  the 
wavelet  network  will  be  too  complex  and  tend  to  overfit  the 
training  data.  On  the  other  hand,  a  wavelet  network  that  is  not 
sufficiently  complex  (over  small  Nw )  can  fail  to  detect  fully  the 
signal  in  a  complicated  data  set,  leading  to  underfitting  [16,24], 
For  this  reason,  the  Cross-Validation  method  is  employed  to 
estimate  Vw .  The  available  data  are  divided  into  two  parts  when 
training  the  model.  One  part  is  used  for  estimating  the  model 
parameters.  The  quality  of  performances  of  the  model  on  the 
other  part  of  the  data  reflects  how  well  it  would  perform  in  an 
unsupervised  setting.  The  Cross-Validation  criterion  is  given  by 

1  N[ 

CV  =  w<*(*))  -  >(*))'  -  (49) 

t  k=l 

where  fyv  is  the  trained  wavelet  network  with  Aw  multi¬ 
dimensional  wavelets.  The  data  set  {(x(k),  y{k))\k  =  1,  ....  N't) 
used  in  Eq.  (49)  is  different  from  the  training  data  set.  Hence, 
/Vw  is  determined  by  minimizing  CV. 

The  modeling  procedures  of  wavelet  networks  can  be  sum¬ 
marized  as  a  flowchart  in  Fig.  3. 

In  the  flowchart,  Aw  should  be  repeatedly  selected  until 
CV  <  ecv-  Several  selection  approaches  exit.  The  immune  algo- 


Fig.  3.  Flowchart  of  the  wavelet  network  modeling  approach. 

rithm  is  a  heuristically  random  searching  algorithm  which  is 
developed  as  an  imitation  of  the  adaptive  humoral  immune 
response  [25],  It  has  been  proved  to  be  capable  of  performing 
such  tough  tasks  as  machine  learning  and  quadratic  optimiza¬ 
tion.  In  this  paper,  we  use  immune  algorithm  for  selecting  Nw  so 
that  the  Cross-Validation  criterion  is  gradually  decreased  after 
selections. 

4.  Results 

Many  factors  such  as  flow  rates,  temperature,  pressure,  mate¬ 
rial  and  configuration  can  influence  the  operating  performance 
of  the  DIR-SOFC.  Furthermore,  the  nonlinear  coupling  vari¬ 
ables  and  parameters  can  increase  difficulty  in  analyzing  the 
operating  states  of  the  DIR-SOFC.  For  a  black-box  model,  we 
need  not  consider  all  the  factors,  and  can  concentrate  on  the 
important  variables  and  parameters  which  we  are  interested  in. 
The  current  density-voltage  characteristic  is  important  for  fuel 
cells;  and  the  thermodynamic  performance  of  the  stack  is  also 
crucial  for  designing  and  controlling  the  DIR-SOFC. 

The  DIR-SOFC  can  be  regarded  as  a  system  with  three-input 
(Fguel,  Fq1  and  V0ut)  and  two-output  (rout  and  /m).  Vout,  Tout 
and  7m  are  the  stack  output  voltage,  the  outlet  temperature  and 
the  mean  current  density,  respectively.  We  assume  that  the  DIR- 
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SOFC  can  be  modeled  by  equations  of  the  following  form. 

T0 „t(*)  =  /i(Vout(*  -  1),  Voul(k  -  2),  F(k  -  1),  F(k  -  2), 

Tomik  -  1),  Tout(k  -  2)),  (50) 

Im(k)  =  f2(Vont(k  -  1),  Vou[(k  -  2),  F(k  -  1),  F(k  -  2), 

Im(k  -  1),  Uk  -  2)),  (51) 

where /i  and/2  are  unknown  nonlinear  functions,  and  F  repre¬ 
sents  ( T0fuel ,  F™)  .  We  expect  to  estimate  functions  f\  and/2  so 
that  the  behavior  of  model  is  close  to  that  of  the  real  DIR-SOFC 
stack.  In  order  to  fit  into  the  framework  of  the  wavelet  network 
as  formulated  in  Section  3,  let 

{x  =  (VoutOfc— 1),  Vou  ,(*  -  2),  F(k  -  1),  F(k  -  2),  Tout  (A:  -  1), 
Tom(k  -  2))t,  y  =  Tout(k)\  (52) 

and 

{x  =  (Voat(k  -  1),  Vout(k  -  2),  F(k  -  1),  F(k  ~  2),  Im(k  -  1), 
Im(k  -  2))t,  y  =  Uk)}.  (53) 

In  our  study,  the  mother  wavelet  has  been  chosen  as 


and  0  as  0.2.  The  wavelet  network  modeling  of  the  DIR-SOFC 
is  implemented  based  on  simulation  data. 

4.1.  Preparing  simulation  data 

The  system  under  consideration  is  a  16  kW  DIR-SOFC  stack 
that  consists  of  30  cells  operating  in  cross-flow  mode.  The 
dynamic  physical  model  of  the  DIR-SOFC  stack  is  developed 


in  MATLAB  to  generate  simulation  data.  The  schematic  illus¬ 
tration  of  the  physical  model  of  one  cell  is  shown  in  Fig.  4. 
The  cell  is  divided  into  several  elements  which  are  mutually 
connected  by  bidirectional  arrows.  The  fuel  and  air  are  respec¬ 
tively  injected  into  the  cell  in  orthogonal  directions  and  pass 
sequentially  through  one  element  to  the  next  element.  The  inter¬ 
action  (temperature,  gas  flows,  etc.)  between  adjacent  elements 
is  taken  into  account  and  is  indicated  by  bidirectional  arrows  in 
Fig.  4. 

Fig.  5  illustrates  the  implementation  of  the  element  in 
MATLAB/Simulink.  The  element  module  is  composed  of  the 
reforming  module,  electrochemical  module  and  temperature 
module.  We  write  S -functions  in  C  language  and  generate 
masked  S-function  blocks  that  can  be  invoked  as  sub-modules.  In 


Fig.  5.  Implementation  of  the  element  in  MATLAB/Simulink. 
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Table  1 

Operating  conditions  and  parameters  of  the  DIR-SOFC  stack 


Item 

Value 

Unit 

Number  of  cell 

30 

Cell  area 

0.16 

m2 

Power 

16 

kW 

Mean  current  density 

100-8900 

Am”2 

Output  voltage 

8-28 

V 

Fuel  flow  rate 

0.5-3 

10-3  kgs-1 

Air  flow  rate 

0.7-3 

10_2kgs_1 

Inlet  fuel  temperature 

980-1320 

K 

Inlet  air  temperature 

980-1320 

K 

Inlet  composition  of  fuel 

CH4 

28% 

moVmol 

C02 

1% 

moVmol 

CO 

1% 

mol/mol 

h2 

10% 

mol/mol 

h2o 

60% 

moVmol 

Inlet  composition  of  air 

o2 

23.7% 

moVmol 

n2 

76.3% 

moVmol 

Operating  pressure 

3 

bar 

Anodic  activation  energy 

1.2  x  105 

Jmol-1 

Cathodic  activation  energy 

1.2  x  105 

Jmol-1 

ym 

2.9  x  108 

Am-2 

yc  a 

7x  108 

Am-2 

Limiting  current  density 

9000 

Am"2 

Fig.  5,  T_ r_out  and  T-neighbour  are  vectors  whose  components 
are  temperatures  of  the  cell  unit,  fuel,  air  and  separators.  T_\n 
and  r_out  denote  the  inlet  and  outlet  temperature  vectors  of  the 
element,  respectively.  Four  T  neighbour  represent  the  temperature 
information  of  the  adjacent  elements. 

The  variables  and  parameters  of  the  dynamic  physical  model 
are  listed  in  Table  1.  The  dynamic  physical  model  is  excited 
with  the  fuel  flow  rate  ((0.5-3)  x  10-3  kg  s-1),  the  air  flow  rate 
((0.7-3)  x  10-2kgs-1),  inlet  temperature  (980-1320  K)  and 
the  output  voltage  (8-28  V)  to  generate  data.  16,000  data  are 
obtained  from  the  simulation  and  are  split  into  two  parts.  One 
part  (8000  data)  is  used  for  the  construction  and  identification 
of  the  wavelet  networks,  and  the  other  part  (8000  data)  is  used 
for  the  evaluation  of  the  resulting  wavelet  networks. 

4.2.  Predicting  with  wavelet  networks 

Using  the  approaches  described  in  Section  3,  wavelet  net¬ 
works  are  thus  constructed  based  on  the  first  part  of  data. 

The  identification  results  of  the  numbers  of  multi-dimensional  2 

0) 

wavelets  used  in  the  networks,  Nw,  are  32  for/i  and  35  for/2  |- 
(corresponding  to  ecv  =  0.3),  respectively. 

Fig.  6  shows  the  comparison  of  V0ut-/m  characteristics  gen-  4= 
erated  by  the  resulting  wavelet  networks  and  the  physical  model.  0 
We  can  see  that  the  static  Vout_fm  characteristics  predicted  by 
the  wavelet  networks  show  good  consistency  with  the  physical 
model  under  various  temperature. 

For  controlling  the  DIR-SOFC  stack,  the  dynamic  predictive 
capability  of  the  model  is  more  important.  Two  cases  are  used  as 
examples  to  evaluate  the  dynamic  performance  of  the  obtained 


Fig.  6.  Comparison  of  Vout-fm  characteristics. 


wavelet  networks. 

(1)  Fuel  flow  rate  increases  from  8.21  x  10-4  to 

9.43  x  10-4 kgs-1,  and  air  flow  rate  maintains 

1.24  x  10-2kgs-1; 

(2)  air  flow  rate  decreases  from  1.24  xlO-2  to 

1.12  x  10-2 kgs-1,  and  fuel  flow  rate  maintains 

8.21  x  10-4kgs-1. 

In  these  two  cases,  the  inlet  temperature  (1024  K),  the  stack 
output  voltage  (21 V)  and  operating  pressure  (3  bar)  are  kept 
constant. 

The  comparison  of  dynamic  responses  of  the  wavelet 
networks  and  the  physical  model  is  shown  in  Figs.  7-10. 
Figs.  7  and  8  present  the  responses  of  outlet  temperature  and 
mean  current  density  for  Case  1,  Figs.  9  and  10  for  Case  2. 

In  Case  1,  the  magnitude  of  changes  in  outlet  temperature 
and  mean  current  density  are  relatively  large,  and  the  curves 
generated  by  the  wavelet  networks  and  the  physical  model  are 
too  close  to  distinguish.  For  a  closer  look  of  the  results,  zooms 


Mean  current  density  (A  rrf 2)  Outlet  temperature  (K)  Mean  current  density  (Am" 
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Fig.  8.  Dynamic  response  of  mean  current  density  for  Case  1. 


Fig.  9.  Dynamic  response  of  outlet  temperature  for  Case  2. 


t(s)  x 104 

Fig.  10.  Dynamic  response  of  mean  current  density  for  Case  2, 


Table  2 

Statistical  results  for  two  i 


Case 

Outlet  temperature 

Mean  currer 

it  density 

MSE  CC 

MSE 

CC 

1 

2 

0.1212  0.9997 

0.0054  0.9992 

0.2106 

0.1189 

0.9993 

0.9989 

Table  3 

Computational  time  for  two  cases 

Case 

Computational  time  (s) 

Outlet  temperature 

Mean  current  density 

Physical  model  Wavelet 

network 

Physical  model  Wavelet 

network 

1 

2 

3.39975053  x  105  0.916 
3.23977547  x  10s  0.872 

3.39975053  x 
3.23977547  x 

1  | 

of  the  curves  are  inset  in  corresponding  figures  (Figs.  7  and  8). 
The  mean  squared  errors  (MSE)  and  the  correlation  coefficients 
(CC)  between  the  wavelet  network  output  and  the  physical  model 
output  in  two  cases  are  computed  and  presented  in  Table  2.  As 
can  be  seen  from  Figs.  7-10  and  Table  2,  the  proposed  method 
is  successful  in  identifying  the  DIR-SOFC  dynamics. 

5.  Conclusions 

Nonlinear  identification  of  a  DIR-SOFC  stack  using  wavelet 
networks  is  described  in  this  paper.  Combining  the  MRA 
of  wavelets  and  evolution  of  networks,  the  wavelet  network 
approach  is  classified  as  black-box  modeling,  and  avoids  using 
complicated  multiple  partial  differential  equations  of  energy 
balance  and  mass  conservation.  The  simulation  programs  oper¬ 
ate  on  a  Pentium  4  2.80  GHz  computer  with  256  MB  memory. 
Table  3  shows  the  comparison  of  computational  time  (required 
for  the  response  curve  to  reach  a  steady-state)  between  the  phys¬ 
ical  model  and  wavelet  networks  for  different  cases  (increasing 
the  fuel  flow  rate  and  decreasing  the  air  flow  rate).  As  can  be 
seen  from  Table  3,  the  computational  time  for  calculating  the 
response  curve  by  using  wavelet  networks  is  less  than  1  s,  and  is 
much  shorter  than  the  computational  time  needed  by  the  physical 
model.  Simulation  results  illustrate  fast  and  robust  adaptation  of 
the  wavelet  network  model  following  changes  in  operating  con¬ 
ditions.  The  static  and  dynamic  characteristics  of  the  DIR-SOFC 
stack  can  be  predicted  with  relatively  high  accuracy.  It  is  shown 
that  the  presented  wavelet  network  is  a  powerful  and  attractive 
method  for  the  identification  of  a  DIR-SOFC  stack.  The  obtained 
wavelet  network  model  can  facilitate  analysis  of  system  stability 
and  be  used  for  designing  model-based  controllers.  In  the  future, 
the  model  is  planned  to  be  used  for  developing  predictive  and 
robust  controllers  for  the  DIR-SOFC  stack. 
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Appendix  A.  Derivation  of  iterative  updating  formulas 

Assume  that,  at  iteration  k,  the  estimation  of  wk  corre¬ 
sponding  to  &is(x(k))  is  u>/s,  while  the  coefficients  of  other 
multi-dimensional  wavelets  are  vj/(k  —  1)  (/  ^  ls).  We  define 

AW 

y(k )  =  Y^wi(k  -  1  )'f,i(x(k))  +  wh  (x(k)).  (A.  1) 

¥h 

A  least  squares  cost  function  with  forgetting  factor  is  defined 
as 

k 

m  =  J2(yU)  -  yU))2^k~j,  (A.2) 


where  0<k  <  1  is  the  forgetting  factor  [23].  wk (k)  is  estimated 
by  minimizing  e(k)  (Eq.  (A.2)).  The  partial  derivative  of  the  cost 
function  with  respect  to  wk  is 


1? 

=  -2j2(y(j)  -  yUms(Mmk~j 

j=  1 

(A.3) 

k  /  Nw 

=  -  2^2  y(J)-^2mU  -  i 

j= 1  V  ¥h 

^(x(j))Xk-j 

-wk&k(x(j))j 

Let 

3e(*)  Q 

dwis 

(A.4) 

Hence, 

wis(k)  =  wis 

£  ■= i  (yU)-Eg$wi(j-  W(x(j)))  *iMmk-j 

Ek=^t(x(mk-j 

(A. 5) 

Let 

dh(k)  =  YJ'l'lCx(j))Xk-J. 

7=1 

(A.6) 

Rewriting  Eqs.  (A. 5)  and  (A.6)  as  recursions: 

A* 

e(k)  =  y(k)  -  ^ wi(k  -  1  )<£/(*(&)), 

(A.7) 

dh{k)  =  kdk(k  -  1)  +  tf(rtk)), 

(A.  8) 

...  n  s(A(k))e(k) 

wh(k)  =  wis(k  1)+  8 

dis(k ) 

(A.9) 

Eqs.  (A.7)-(A.9)  are  used  for  updating  w\r 
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